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1. Introduction 

Condensed states of identical bosons at zero temperature in a confining external field 
are accessible for experiments in many laboratories Qllll- Most often the condensates 
are formed by one species of bosons, but also two-component structures have been 
investigated. The two components may consist of the same atom in two different spin 
states as studied for ^"^Rb jSl II]. The spin states may then be converted from one to 
the other while maintaining the total number of atoms. Also collective oscillations can 
be studied experimentally as in [S] for ^^Rb. Two-component boson condensates of 
different atoms have other properties as investigated by combining ^^Rb and ^^K [Hj. 
More recently mixtures of boson and fermion systems have been realized [3 IHl 13 ^1 ^] • 

The theory for two-component condensates are mostly based on the mean-field 
approximation ^ The additional degrees of freedom lead to spatial symmetry 
breaking in a cylindrical trap as shown in ^21 E]- Collective excitations in such traps are 
predicted ^1]. With two different trap frequencies also spherical traps allow collective 
excitations ^S]. The spatial asymmetry is sometimes favored even in spherical traps 
Repulsion between the two components can lead to separation of the two components 
both in the ground state and for the collective oscillations jTJj . Stability conditions can 
be derived from properties of the collective modes [TH] . 

Another type of formulation using hyperspherical coordinates for one-component 
boson systems was introduced in This formulation was extended to include two- 
body correlations for one-component iV-body boson systems This treatment 
allowed both very large scattering lengths and attractive two-body interactions. The 
mathematical collapse is prevented by a finite range interaction. Conditions for the 
physical collapse are rather similar to mean-field results. Still the computed correlated 
structures for large scattering length may be physically interesting pTj . 

The present paper follows up with more details and new calculations on a recent 
attempt to extend the hyperspherical formulation to two-component boson systems 
[22] . Although several new features appear, we still only employ the lowest-order 
approximation with inclusion of three collective coordinates. This is the minimum 
number of degrees of freedom which can describe individual sizes of the two subsystems 
and their separation distance. A demand for more complicated extensions may arise, 
but the specific direction is better decided after digestion of the results from the simplest 
model survey. 

In section 121 we formulate the theoretical framework where the effective potential is 
introduced. The properties of this key quantity is investigated in details in section 
El especially with respect to minima and their curvatures. Then we are equipped 
to estabhsh stability conditions for the total system as discussed in section HI Our 
conclusions are contained in sectional Finally, for completeness we include a number of 
rather tedious mathematical details in three appendices, but the main text can be read 
independently because only analytic and numerical results are used. 
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2. Theoretical formulation 

We consider a dilute system of two species of bosons, A and B, with boson masses ttia 
and mB and particle numbers A^^ and Nb- We label the particles with i = 1,2, 3, Na 
for compontent A and i = Na + 1, Na + 2, ...,N = Na + A''^ for component B. The 
individual masses and position vectors are denoted and where i = 1,2,3, N. 
These bosons are trapped in a spherically symmetric harmonic oscillator potential with 
the same trap frequency u. The Hamiltonian describing this system is 

H = Ha + VA + VB + Vab (1) 
- ^ + t l-^V? + l,nB^r?] (2) 

1=1 i=NA+l 
Na . »-2 

Va= E^^^^'Hn.), 9A = ^^, (3) 

, TTlA 

i>j=l 

Vb= 9BS^''(r^,), gB = ^-^, (4) 

Na N ,2 

Vab = Y: E 9ABS^'\r^,),gAB=^-^^^, (5) 

where Vij = r j — Vj denotes the vector distance between two particles and the reduced 
mass of an A and a B particle is given by fj,AB = t^AiT^Bl ij^A + ^b)- The interaction 
strengths of the two-body zero-range pseudo potentials, qa, Qb, and qab, are related to 
the s-wave scattering lengths, oa, a^, and cab- This normalization is adopted to provide 
the correct large- distance behavior of the effective potential in mean-field computations 

We shall use hyperspherical coordinates with hyperradii defined for both the 
individual and the total systems [1111201, i-e. 

^ N N 

i<j i=l 
^ Na 1 ^ 

Kj NA<i<j 
^ Na 1 ^ 

R^ = Wa^''^ ^' = Wb S 

^ i=l ^ i=NA+l 

where M = Ma + Mb = Nattia + Nb^b is the total mass, m is an arbitrary 
normalization mass, Ra-, Rb are the individual center of mass coordinates and R is 
the overall center-of-mass coordinate, i.e. 

R=LYm-r. = ^l^^A±^l^ (9) 
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The center of mass separation of the two components is then given by the coordinate 
r = Rb — Ra- The three length coordinates, pa, Pb, and r, are related to p by 

2 2 , 2,2 _ MaMb , . 

mp = tuaPa + ^bPb + , mr = — — — . (10) 

Ma + Mb 

The remaining coordinates, beside pA, Pb, R, r, are all chosen as (hyper) angles denoted 
collectively for each system by Qa and Qb- 

The total volume element is then given by 

Na N 

i=l i=NA + l 

d^Rr'^drdpAdpBdflrdf^^NA-id^NB-i y (H) 
where dflr is the angular volume element for r and dflj\f^_i,dfl]\fg_i accounts for the 
hyperangles collected in Qa and Qb, see [201 for precise definitions. 

Using these coordinates we rewrite the external harmonic oscillator potentials in 
the Hamiltonian as 

-niAUJ^ ^ri + -niBUJ^ ^ = 

1=1 i=NA + l 

^Muj^R^ + im^cjV^ + ^TjiAUJ^pl + ^niBUJ^pl (12) 
2 2 ^ 

Furthermore the kinetic energy operator T can also be separated in terms related to 
total center of mass, relative and intrinsic A and B degrees of freedom, i.e. 

where V_r and Vr- are the usual three-dimensional derivatives with respect to the 
indicated coordinates, Ta and Tb are the intrinsic kinetic-energy operators expressed 
by hyperspherical coordinates for each boson system, see [20]. 

The two-body interactions depend only on relative distances and we can completely 
separate relative and total center of mass motions. Furthermore, we can separate relative 
A and B except for the coupling term Vab- The Hamiltonian in equation then 
becomes 

H = Hem. + Hj-ei , (14) 

Hem. = -^fl + \Muj'R\ (15) 

111 
i^rei = Tp^ + Tp^ +Tr + -niAUJ^ p\ + -^^^^^1 + -m.cuV^ + H^^ (16) 

B^ = ^ + ^ + ^ + Va + Vb + Vab , 17 

2mA Pa 2m b Pb 2mr 

1 92 

Tr = ' (18) 

2mr r or'' 
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T - J-3Na/2 3Na/2-2 /.qx 

T - 37Vfl/2-2 /^^v 

2mB^^ dp'/"" ' ^^"^ 

where is the usual angular momentum operator with respect to r, and A^^_]^ and 
consist of first and second order derivatives with respect to the hyperangular 
degrees of freedom indicated by the indices. This Hamiltonian reduces for ttia = itib 
and identical interactions, qa = 93 = Qab, to the one-component system described in 



The separable center of mass motion is given by the harmonic oscillator solutions 
to the Hamiltonian in equation pSj). The relative motion can be determined by the 
hyperspherical adiabatic expansion method where the wave function for fixed values of 
{PA,PB,r) is expanded on the complete set of eigenf unctions for Hq in equation (fT7|) . 
The leading term in the lowest adiabatic channel is expected to consist of the lowest 
partial waves, where the wave function is independent of all (hyper) angles. This is the 
assumption valid in the large- distance limit where all particles are far from each other 
and all directional dependence is averaged out [IHl El! ■ Then the wave function for the 
lowest eigenvalue reduces to the form 

,j, 2-3Na/2 2-3Nb/2 _W/ x \ 

^ = Pa Pb ^ f{pA,pB,r), (21) 

where / only depends on the three radial coordinates. 

This radial wave function and the corresponding eigenvalue i^rei are determined 
from the Schrodinger equation obtained by combining equations and ()2H] . i.e. 



2mA dp\ 2m b dp\ 2m.,. dr'^ 

+ f/efr(pA, Pb, r)]f{pA, Pb, r) = E.^ifipA, Pb, r) (22) 

UcsipA,pB,r) = -u'^{mAp\ + -rriBpl + -my^ 

n^(3A^A-4)(3A^A-6) ;i^(3iVB -4)(3iVB -6) 
8mAPA 8mBP% 
+ {%\Va\%) + {^o\Vb\^o) + {^o\Vab\%) , (23) 

where ($o|^|'^'o) is the expectation value of the interaction V with the constant angular 
wave function $o- With the volume element in equation (jlip we find by direct integration 

($o|l^A|«fo) = -^7T71^Na{Na-1)—4, (24) 



X ^Na{Na - msiNa - 1)I(PA, PB, r) , (26) 
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where T{z) is the Gamma function and I is derived in Appendix A as a function of p^, 
Pb and r. 

When iV^ > 1 we have r(^^^)/r(^^^) ^ (3iV^/2)3/2 analogously for 
A'"^ ^ 1, see . Then Urmeff is given as 



1 2 2 1 2 2 1 2 2 mlh^ 

U^, = -m^up, + -ms.p, + -m,ur + ^ + ^^^^ 

33/2;^2 , Ar7/2^ . atV^ , 



+ 



which can be considered as the energy surface depending on the three collective 
coordinates p^, Pb and r. The properties of this effective potential then reflects the 
properties of the solutions which, if necessary, could be computed from the eigenvalue 
equation (j^ . 

The expression reduces to the one-component potential [H EDI when r = and all 
particles are identical, i.e. 

1,2 QAT^n^ ?,^l^h^Ny^as , , 

t/cff = -mu p + Y + T^TTT^ — ' 28 

2 8mp^ 2"'7r^/"^ mp'' 

where is the common scattering length and m = tua = ^b- 



3. The effective potential 

To be specific we consider parameters corresponding to a system where both components 
are '^^Rb atoms with the mass itla = niB = itl = 1.44 x 10~^^kg The bosons in both 
components are assumed to be trapped in the same spherical harmonic potential with 
the frequency = 27r x 23.5Hz. The related oscillator length bt = y^h/muj is then 2.22 
pm. The particle numbers N^, Nb and the three s-wave scattering lengths a a, as and 
ttAB are variable parameters. 



3.1. The interaction Vab 

All terms in the potential ()27|) are simple and explicitly given, except the function / 
arising from the interaction Vab in equation p6|l , see appendix [Appendix A[ Without 
this interaction the two subsystems would decouple and each behaves as a separate one- 
component system. We show this interaction term in figures ^ and El for several sets of 
particle numbers and hyperradii. 

The general features are a monotonous function with a fiat maximum at r = and 
vanishing exactly for r > pa + Pb- For given A^ the symmetric combination of Na = Nb 
and Pa = Pb leads to the highest maximum value of unity corresponding to the smallest 
average distance between particles in the subsystems A and B. The interaction between 
the two components is strongest for symmetric systems with coinciding centers. As the 
centers are moved apart, the interaction becomes larger for different values of pa and p^. 
Asymmetric division of a given total number of particles also decrease the interaction. 
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Figure 1. The interaction / defined in equation 1)26(1 and [Appendix A| as function 
of the distance r between the two centers of masses. The particle numbers are 
Na = Nb = 20000, Pa = 1006t and the ratios pb/pa are given on the figure. 



These properties directly reflect the dependences of the overlap between the two density 
distributions. 
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Figure 2. The interaction function / defined in equation (|26|l and [Appendix A| as 
function of the distance r between the two centers of masses ioT pA = Pb = 1006t and 
the values of {Na, Nb) given on the figure. 



The two-body interaction strength qab is a proportionality factor on the function 
/. Consequently, attraction prefers symmetric division and coinciding centers, whereas 
repulsion prefers precisely opposite, i.e. asymmetry and large distance between the 
center of masses. 

3.2. Dependence on center of mass distance 

We minimize the effective potential in equation (|7fj) with respect to pA and ps for 
fixed r and given scattering lengths and particle numbers. We first choose interactions 
corresponding to the two-component condensate of two spin states of ^''Rb. The three 
scattering lengths are then almost identical and given hj qa = clb = clab = lOSag where 
ao is the Bohr radius The total number of particles is rather arbitrarily fixed to be 
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N = Na + Nb = 40000. The results for identical scattering lengths are discussed in 

The main features are that the effective potential at r = is independent of the 
division of particles and Nb, but still strongly varying with the total number A^. 
This potential value at r = is the same as obtained for the one-component potential 
in equation (^Hjl when p is calculated from equation (fTUIl with the minimum values of 
Pa and pb- 

For these identical repulsive scattering lengths and the symmetric divisions of 
A'"^ = A^^, the systems prefer coinciding centers of mass, i.e. the potential has a 
minimum for r = 0. As we increase the particle asymmetry, another minimum appears 
at r ^ 2 — 3bt and for Nb < 0.134A^^ (or A^^ < 0.134Nb) eventually becomes deeper. 
The minimum at r = remains, but now separated by a barrier. The deepest of these 
minima is found when Nb/Na ~ 1/40. Even larger asymmetries again increase the 
depth and eventually only the r = minimum remains as for a one-component system. 

Maintaining = as > we show in figure El the variation with the scattering 
length qab > for a symmetric division of particle numbers. The flat region around 
r = is always present, but the minimum for aAB < cla turns into a maximum when 
ttAB > O'A- Then the repulsion between the subsystems is strong enough to move the 
minimum at r = to finite values corresponding to values r/bt ~ 3.5 — 4.5. For less 
repulsion, aAs < o-A = clbi the minimum at r = is always present and in addition 
another minimum at finite r appears for sufficiently asymmetric particle divisions. 




1 2 3 4 5 6 



Figure 3. The effective potential in equation (|17ll as a function of r/ht for ^^Rb masses 
■niA — ruB = 1.44 x lO^^^kg, trap length bt = 2.22/im = 27r x 23.5 Hz), particle 
numbers iV^ = Nb = 20000, and ua ^ clb = 103ao, where uq is the Bohr radius. The 
curves correspond to the different values of gab I o-a given on the figure. 



The basic characteristics of a system are size and energy. The root-mean-square 
distance is the measure of the size of a system which in our coordinates are given by 
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i=NA + l 

N 



AT, 



B 



{n - Rf) 



M 



(30) 



(31) 



where f, ta and tb are the root mean square radii for total and individual subsystems 
and represents the expectation value for the total wave function. 

For identical interactions of qab = cla = as found in j2Z| that = = r ^ 3bt 
in the minimum at r = independent of particle division. This structure closely 
resembles the one-component system. Increasing the distance r between the centers 
of the subsystems reduce all three root mean square radii. However, as Nb/Na changes, 
the system with the largest particle number always maintains roughly the same value 
of about 3bt. In contrast, the root mean square radius of the other subsystem decreases 
to about half of the value for r = 0. 

The minima of effective potentials as in figure El define the preferred structures. The 
effect of the repulsion is that the two systems try to avoid each other while staying at 
distances small compared to the size of the confining external field. The large subsystem 
is exploiting the space almost like it was alone in the trap. When one subsystem is small, 
it becomes advantageous to place about half of its particles outside the other subsystem. 
Then the repulsion on the small subsystem from the trap and the other subsystem is 
minimized. 
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Figure 4. The root-mean-square radii ta = as a function of r for some of the 
systems in figure 01 The curves again correspond to the different values of a^s/a^ 
given on the figure. 



When Nb < 0.134A^a, the lowest minimum occurs at finite r. The center of mass of 
the small system B remains within the radius until Nb ~ O.OOANa- For even smaller 
subsystems the center of B moves a little outside f^, but still the distance between the 
centers at the minimum remains smaller than fA + tb- Therefore, complete separation 
between the systems does not occur. 

Many of these basic properties remain for < qab cla = cib &s shown in figure 
m The root mean square radii again are largest for r = with values increasing as a 



Two- component boson systems with hyperspherical coordinates 



10 



function of a^B corresponding to decreasing the overlap in order to decrease the energy. 
For relatively large values of r where the subsystems essentially can avoid overlapping, 
the sizes become independent of aAB- At intermediate distances a minimum size is 
always found. The smaller the repulsion a as, the less is the variation with center of 
mass distance. 



3.3. Curvature and normal modes 



The local minima in the r-direction represent a projection in the three dimensional 
space. We can evaluate the stability by computing the second derivatives of the effective 
potential in these points. The kinetic energy operators already only contain second 
derivatives and the relative Hamiltonian H^ei is then to second order approximated by 



H^cl ~ UcsipAmin, PB 

2mA dp\ 
1 

+ 2 



mm; ' mm ) 



2mB dp% 



2mj. dr"^ 



Pa — Pad 



( d^u, 

— — Pi 



PB ~ PBmm, 



r — rr 



dpBdpA 



dpAdpB 



y drdpA 



9^ 



drdpB 



dpAdr 



dpB dr 



( Da — 



J 



\ 



Pa ~^ PAm'm 
PB — PBmin 
^ ^min 



(32) 



where we used the notation of the 3 x 3-matrix sandwiched between the two vectors. 
All first order derivatives are zero in these minima. 

The approximation in equation (jH^ is almost directly three separable harmonic 
oscillators. We first redefine the variables by scaling with the related masses as 

^PA = {PA - PAm\n)\/mAi^/fi , (33) 
SpB = {PB - PBmm)\/mBU^/h , (34) 

6r = (r - rmin)A/r?v^ • (35) 

Then all three kinetic energy operators are dimensionless and the energy unit is hv. 
It is now sufficient to diagonalize the potential matrix in equation ()32j) . The total 
Hamiltonian then has eigenvalues and eigenvectors Aj and vi = {vf, vf, f[) for i = 1,2, 3. 
The total energy is given hj E = hwY^^i^i + 1/2) with the non-negative integers 
Ui corresponding to one-dimensional harmonic oscillators. The excitation energies for 
each of these vibrational degrees of freedom are then integer multiples of ^a/A^. 

For the minimum at r = the pA-pB and the r degrees of freedom decouple 
completely. This is a consequence of the flat maximum at r = for the interaction 
function in figures Q and |21 because first order r-derivatives of / then vanish and the 
only other r-dependence is the arising from the external field. For this minimum the 
lowest vibrational r-mode carries the energy fijuj{nr + 3/2)y/X^ dictated by the boundary 
condition at r = corresponding to negative parity. The total energy is in this case 
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Table 1. Eigenvalues and eigenvectors for the second order harmonic motion in 
the minimum at r = of the effective potential. The eigenvectors {vfjvfvD refer 
to the units of root mean square radii or in terms of the initial coordinates {{pA — 
PAmin)\/NAmAUj/h, (fB - R Bmin) \/ N sm B^J / K " r^in) \/ mrUj / H) . The scattering 
lengths are uab — cla — o-b — lOSao- 



Na 


Nb 


Ai 






vl 


A2 


V2 


v^ 


vl 


As 






vl 


20000 


20000 


4.97 


0.71 


0.71 


0.0 


2.54 


0.71 


-0.71 


0.0 


0.03 


0.0 


0.0 


1.0 


32000 


8000 


4.97 


0.97 


0.24 


0.0 


2.54 


0.71 


-0.71 


0.0 


0.03 


0.0 


0.0 


1.0 


36000 


4000 


4.97 


0.99 


0.11 


0.0 


2.55 


0.71 


-0.71 


0.0 


0.03 


0.0 


0.0 


1.0 


37000 


3000 


4.97 


1.00 


0.08 


0.0 


2.55 


0.71 


-0.71 


0.0 


0.03 


0.0 


0.0 


1.0 


39000 


1000 


4.97 


1.00 


0.03 


0.0 


2.55 


0.71 


-0.71 


0.0 


0.03 


0.0 


0.0 


1.0 


39800 


200 


4.97 


1.00 


0.01 


0.0 


2.56 


0.71 


-0.71 


0.0 


0.04 


0.0 


0.0 


1.0 



given hy E = fvjj{{nz + 3/2)\/A^-|- + l/2)"\/Ai)- Still the vibrational excitation 

energies are integer multiples of tvjjy/Xi. 

The normal modes are given by the components of the eigenvectors Vi defining 
a direction in the coordinate system {6pA,SpB,Sr). The physical system has natural 
length scales for the different degrees of freedom corresponding to the initial coordinates 
{PA,PB,r)- They are exhibited by the sizes of the subsystems, i.e. the root mean 
square radii defined in equations (f^ - (jHT|) . Thus the proper physical length scales 
are pa/V^a and pb/V^b, which means that we should multiply the coefficient 
measuring the p^-component by \/Na to change to the coefficient measuring the 
corresponding root mean square radius. The normal mode directions with these 
root mean square unit vectors are then represented by the vectors {vf, vf, v^) = 
{vf ^/NAmAUJ/h, vf y/NBniBu/h, vl ^Jrn^wjK) . 

These eigenvectors and the related eigenvalues are given in tables Q and |21 For the 
minimum at r = the r-direction is completely decoupled from the two other normal 
modes mixing pA and pb- By far the lowest eigenvalue corresponds to the r-mode which 
means vibrations of the relative center of mass. 

The next mode is remarkably independent of the particle division, always with 
exactly equal amplitudes of opposite phase in the A and B root mean square radii. This 
mode is then an isovector breathing mode where one subsystem shrinks radially while 
the other radially expands. This vibration is around the equilibrium structure where 
the subsystems both are of equal size. 

The highest-lying mode for the r = minimum corresponds to in-phase oscillations 
of Pa and pb- The equal amplitudes of the real size components, (1,-1) oc 
{v^-,^^^ oc (f^A/iVyi, V2\/Nb), combined with orthogonality of (f^,ff) and (^2^,^^) oc 
{y/Ns, -VNa), implies that {v^,v^) oc (v^, v^), and {v^,v^) oc {Na, Nb). 

In general, the vector {v^, v^, v"") oc (pAminxA^jTW^, PBm:m\frnB^jJh, r^^^^uirUJ /h) 
is at the minimum position precisely in the direction of the total hyperradius defined 
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Table 2. Eigenvalues and eigenvectors for the second order harmonic motion in 
the minimum at r 7^ of the effective potential. The eigenvectors {vfjvfvD refer 
to the units of root mean square radii or in terms of the initial coordinates {{pA — 
PAmin)\/NAmAUj/h, (fB - R Bmin) \/ N sm B^J / K " r^in) \/ mrUj / H) . The scattering 
lengths are uab — cla — o-b — lOSao- 



Na 


Nb 


Ai 






vl 


A2 


V2 


vi 


vl 


A3 








20000 


20000 


























32000 


8000 


























36000 


4000 


4.97 


0.99 


0.08 


0.08 


3.29 


0.72 


-0.62 


-0.31 


0.26 


0.26 


0.49 


-0.83 


37000 


3000 


4.97 


1.0 


0.06 


0.06 


3.47 


0.71 


-0.62 


-0.34 


0.48 


0.32 


0.51 


-0.80 


39000 


1000 


4.97 


1.0 


0.01 


0.02 


3.81 


0.67 


-0.63 


-0.40 


1.21 


0.45 


0.50 


-0.75 


39800 


200 


4.97 


1.0 


0.0 


0.0 


3.90 


0.60 


-0.67 


-0.44 


1.77 


0.53 


0.47 


-0.71 



in equation (|10p. For r = the root mean square radii of the two subsystems are 
equal and we have from equations (j^ - lj^ that pAminV^s ~ PsmmV^A- Thus the 
highest-lying vibrational mode, (f ^, f f ) oc {^/Na, V^b), is exactly in the direction of 
the hyperradius. 

For an equal particle division this is fully the physical breathing mode where both 
subsystems move in phase with an equal scaling of their radii. For asymmetric particle 
divisions the largest subsystem tries to breathe like it was alone. The smallest then has 
to follow as well as possible consistent with orthogonality, which in turn is equivalent 
to following the p-direction. Equal breathing amplitudes are then only consistent with 
the p-direction for a symmetric particle division. 

For the minimum at r 7^ all three directions are mixed in the normal modes 
as seen in table El The lowest-lying energy is now much larger than for the r = 
minimum, but the largest probability is still related to vibration in the r-direction. The 
second mode again corresponds almost completely to the isovector breathing mode, i.e. 
out-of-phase oscillation with equal real size amplitudes of the two subsystems. Now an 
admixture of up to 30% probability in the r-direction is present. 

The highest-lying mode is in-phase oscillation of the two subsystems. Again this 
mode corresponds to maximum variation of the total hyperradius and therefore the 
breathing mode cannot be fully exploited. The reason is that the largest subsystem 
determines the mode and tries to oscillate as it was alone in the trap. 

For both minima the smallest eigenvalue corresponds essentially to the r-direction. 
This mode exploits that only the repulsion between A and B is changing when the 
subsystems maintain their sizes and only move their centers of mass. The second 
and third eigenvalues are essentially the isovector and isoscalar breathing modes, 
respectively. The highest vibrational excitation energy, hhj^/X, is consistently about 
hLO\f^ which is the one-component result for a Bose-Einstein condensate in the limit 
of large particle numbers |2]. This breathing mode is energetically less favored than 
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the isovector breathing mode, where one subsystem oscillates against the other and the 
repulsive interaction minimizes the overlap. In both these modes all three repulsions 
contribute. 



4. Stability conditions 

The stability is first of all determined by the existence of minima. However, for many- 
body atomic systems lower lying minima certainly exist with structures completely 
different from condensed states. These metastable states are confined by barriers 
providing a finite lifetime. In three-dimensional quantum mechanics the minimum in 
the potential energy may be too shallow or narrow to allow the zero point oscillations. 
Crude estimates of stability are then first obtained by locating extremum points with 
positive curvature in the potential. Second by comparing the vibrational energies with 
the barrier heights. 



4.1. The r-mode for identical interactions 

Let us first assume identical scattering lengths where two minima in the r-direction 
appears for some divisions of particles Na and Nb- The values of the effective potential 
at these two minima for r = and r 7^ are shown in figure El The first prominent 
feature is that only the r = 0-minimum exists for Nb/Na > 5400/34600 ~ 0.156. For 
more asymmetric divisions, Nb/Na < 4712/35288 ~ 0.134, the minimum at finite r 
quickly becomes substantially deeper. 

The zero point energies are added to the potential, but the sum coincides with the 
potential within the thickness of the lines in figure El Both oscillations are essentially 
always classically allowed, i.e. around finite r very quickly after it appeared and around 
r = until the asymmetry reaches very small values. 



3 

% 7.278 - 




H 7.268 -, 



Nb/Na 



Figure 5. The values of the effective potential at the extremum points as a function 
of Nb/Na for Na + Nb — 40000 and qa — cib — o-ab — lOSao- We show the minima 
at r = (solid line) and r 7^ (long-dashed line) , and the separating maximum (short- 
dashed line). The zero point energies, added to the effective potentials, are also shown 
for both minima, r = (dotted line) and r (dot-dashed line). 
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The highest minimum located at r = could in principle decay into the lowest 
minimum at finite r. The related lifetime can be estimated by the WKB tunneling 
approximation where the decay rate T/h is given by 

T/h ='±exp(^--J^ dr^/2mr{Ues{r) - U,s{r = 0) - 3;Lu,/2)) , (36) 

where ri and r2 are the turning points and Ur = UJ^/X^■ is the oscillation frequency in 
the r-direction. These decay rates are computed and shown in figure IHl for the r = 
minimum for the cases exhibited in figure El The lifetime increases dramatically as 
the systems become more asymmetric. Only when Nb/Na is close to 0.001, the rate 
becomes comparable to the frequency of the trap. For example for Nb/Na ~ 0.0025 
the rate T/{nw) ^ 8 x lO^^. 

10^ p~i — I — I — I — I I 1 1 1 1 — I — I — I — I I 1 1 1 1 — I — I — I — I I i_ 




0.01 0.1 1 

Figure 6. The WKB decay rate defined in equation H36fl as function of Nb/Na for 
the minimum at r = for the cases in figure [S] 

The structures in these minima are refiected by the root mean square radii defined 
in equations ()29|) - (jHT|) . For r = all these sizes are equal. For finite r we show the results 
in figure[7|for the same parameters as in figureEl The largest system maintains the same 
size ta ~ 3^4 independent of particle number. The size of the smallest system decreases 
with increasing asymmetry. In all cases we find that ta+tb > f, i-e. complete separation 
between systems never occurs. Furthermore, f ^ > f when Nb/Na > 150/39850 ~ 0.004, 
i.e. the center of mass of the small system remains within the radius of the large system 
except for very asymmetric particle divisions. 

4-2. The r-mode for symmetric particle division 

The stability of the r-mode changes with the three interactions. For r = the condition 
for instability against separation of the two center of masses is a negative curvature 
in the minimum of the effective potential in the r-direction. We restrict ourselves to 
tua = mB = m, Na = Nb and aA = clb and define the coordinate at the potential 
minimum by pAmin = PBmin- This minimum point is determined as the zero point of the 
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Figure 7. The root mean square radii ta (solid), (long-dashed), r (short-dashed) 
and fA + tb (dotted) for the minimum at finite r shown in figure [S| The minimum at 
r — has fA=rB~ 36t which also is shown for the most symmetric systems. 



first derivative, i.e. ^^1^=0 = 0, which is equivalent to 

f{pA) = [muj^p\ - - «A + aAB)] = , 37 

with precisely one solution pAnan/bt > (3iVA/2)i/25-i/4 when iV^(aA + aAB)lht > 
— (87r)^/^5~^/^. Furthermore, the function /(pa) is positive for pA > PAmm- 
A negative curvature, \r=o < 0, is then equivalent to 

< 27rV2 ^2^2 = Pac (38) 

The instability condition, (3A^^/2)^/25"^/^6t < PAmin < Pac, IS equivalent to the 
inequahty f{pAc) > 0, which in turn is equivalent to the two simultaneous conditions 

r^foABY' 
bt ^ h 2y^Nf \ h ) 



aA ^aAB_ , (39) 



> (40) 
bt 23/255/4 • ^^^^ 

When = ctB = the critical repulsion derived from equation (pUj) corresponds 
to aAB/bt = 2^Na ' "^"^^ fiiiite values of = ct^ > the critical value of aAB has to be 
correspondingly increased. The behavior of the effective potential is illustrated in figure 
El for parameter combinations close the critical values for separation of the two centers 
of mass. The potentials are rather fiat and the minimum therefore quickly moves to 
separation distances comparable to the trap length. 

Ji-.3. Stability of the pa and pB-modes 

Qualitative understanding can be reached by analytic investigations for symmetric 
systems, i.e. Na = Nb, aA = aB with the constraint r = 0. Then the condition 
that the effective potential at r = has no stationary point, ^j^\r=o 7^ for any 
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coordinates pa and pb, can be expressed as 
NA{aA + aAB) - 23/2vrV2 



< 



-0.67 . 



(41) 



bt - 55/4 

Then the two-component system is unstable since collapse to a point is unavoidable. 

Stability requires, beside a stationary point, also that the curvature is positive in 
all directions. Then the eigenvalues of the curvature matrix must all be positive at the 
stationary point. For r = we then construct and diagonalize the matrix 

dpsdpA dpg / 

For Pa = Pb can compute the eigenvalues analytically and find the conditions 
expressing when at least one of them is negative and the point therefore unstable. In 
Appendix B| we derive these conditions, i.e. the system is unstable when the following 
two inequalities both are satisfied 

5 



CtAB 

NAaA 



< 



< - 



23/2^1/2 



OA 

bt 



+ 



Ioa 
4 bt 



(43) 
(44) 



The last inequality, the stability condition for the A-component alone, can be 
obtained both in the mean-field approximation and with hyperspherical coordinates 
PI 121 E]- Furthermore, equation (j44p is obtained from equation (j43j) for non-interacting 
subsystems, where qab = 0. Both equations (PT|) and PHj) can be derived in the mean 
field approximation HF 
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Figure 8. Contour plot of the potential Ucfi / {Nfiio) in equation H27() as a function 
of PA and PB for r ^ Na ^ Nb ^ N/2 = 2000, oab = 0.00016*, and 
aA=aB= -0.000426f. 



These conditions are most easily visualized by two-dimensional contour diagrams 
of the effective potential as function of pa and ps with r = 0. Such a diagram is 
shown in figure |H1 for interactions very close to the critical values given in equations 
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(j43|) and (j44j) . The potential increases when either oi ps individually become large. 
This effect of the external field is especially clearly seen along the line pa = Pb, where 
the potential decreases from large values at very large distances, proceeds through a 
minimum located at around 426^ to reach a maxium at about 306^, and finally decreases 
to — oo at the origin. However, the minimum along this line is in fact a saddle point and 
therefore unstable as soon as both degrees of freedom are included. From the saddle 
point the energy decreases in the "isovector" direction where the sum of pa and ps 
remains constant. 

This saddle point changes into a stable minimum with slightly larger repulsion 
between the two subsystems as illustrated in figure IHl The same qualitative behavior 
is seen along the Pa = Pb line, where the minimum and maximum are shifted to 
positions further apart at about ASbt and 206*, respectively. Also the regions at relatively 
large distances along the coordinate axes are essentially unchanged. In the "isovector" 
direction, as well as in all other directions, the potential now increases as characteristic 
for a minimum. However, the minimum is very shallow and only a small amount of 
energy would destabilize the system. 
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Figure 9. Contour plot of the potential UoS / {Nfiio) in equation 1)2 7|l as a function 

of PA and PB for r ^ Q, Na ^ Nb = N/2 = 2000, aAB = 0.0002&t, and 
aA=aB= -0.00042fot. 

4.4- Phase diagram of stability 

We are now in a position to discuss the stability for different interactions. We 
consider again equal particle division while independently varying the scattering lengths 
OA = as < and aAB- The numerically computed stability diagram is shown in figure 
ITUl We shall discuss the different regions in this contour plot in comparison with the 
analytic derivation. 

First, we consider two attractive subsystems, aAB/h < 0, where the coinciding 
centers are preferred. The computed stability condition then follows the expression in 
equation (jUj), which for aAB = reduces to the one-component stability condition. 
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Thus, the attraction between A and B reduces the stabihty region. It is amusing that 
= = leads to the apparently identical stability condition of NAaAs/h > —0.67. 
However, since Na = Nb now the total number of pair interactings is twice as large, i.e. 
N\ compared to 
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Figure 10. Stability regions for a two-component system with Nb = Na = 2000 as 
a function of = < and uab ■ The thick solid line separates stable (right side) 
from unstable (left side) regions when we allow three degrees of freedom, i.e. pA, Pb, 
and r. The different curves are related to equations 141(1 (dot-dashed), 143|l and 1(44(1 
(long-dashed), ((39() and 1(40(1 (short-dashed), respectively. 

For relatively small repulsion, < aAs/h ^ 0.0001, the stability follows the 
expression in equation PHj) . This corresponds to a local, but unstable minimum at 
r = 0, since equation then also is obeyed. These conditions are related to structures 
of coinciding centers and therefore they are the same as derived from the mean-field 
approximation [TS] . 

For larger repulsion, 0.0001 < aAs/bt, separate center of masses of the two systems 
are energetically favored. The condition for separation is given in equation ()40|) which 
follows the computed stability curve in the interval 0.0001 < aAs/h < 0.0003. The 
substantial reduction of the region of mean-field stable solutions is then defined by the 
difference between the curves described by equations (pn|) and (pSjl . 

For larger repulsion aAs/bt > 0.0003 the stable region is larger than found 
from equation ()40|) . because the structure at finite r is stable for relatively strong 
repulsion between the subsystems. This repulsion has the effect of squeezing each of the 
subsystems more than if they had been left alone as described by equation (0^. Then the 
stable region in figure ITUl is in fact smaller than for non- interacting subsystems, because 
the barrier preventing collapse now is circumvented by the energetically favored smaller 
sizes and the finite r. Each subsystem therefore finds itself with values of pAmin and 
PBmin on the unstable side of the barrier preventing collapse if each were left alone. The 
mutual repulsion destabihzes the total system compared to separate subsystems. 
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5. Conclusions 

We use hyperspherical coordinates to describe a mixture of two different components 
each consisting of identical bosons. Sufficiently dilute and spatially extended systems can 
be rather well approximated by use of only s-waves for all relative angular momenta. 
From far apart only the monopole part of the interactions is important and this is 
reflected in the wave function. We are then left with a wave function only depending 
on the most important coordinates, i.e. the hyperradius for each subsystem and the 
distance between the two center of masses. 

With a wave function depending on these three collective coordinates and three- 
dimensional (5-functions as the two-body interactions we derive the corresponding 
Schrodinger equation and compute the effective potential. The behavior of this potential 
and the external field as a function of the three collective coordinates are decisive for the 
properties of the two-component system. Several features are similar to results from the 
mean-field approximation. We used the simplest non-trivial assumptions and leave lots 
of room for improvements. Especially two-body correlations could be included in the 
wave function and finite range two-body interactions with the correct behavior could be 
employed. 

When the interaction between particles in the two subsystems is attractive, the 
potential has a minimum for coinciding centers of mass. Then the stable structures are 
quite similar to those obtained with the mean-field approximation. When the repulsion 
is sufficiently strong, new stable structures arise with finite distance between the centers. 
Asymmetric particle numbers in the subsystems favor these structures. For moderate 
repulsion two stable structures are simultaneously present for centers coinciding and at 
a finite distance. 

New vibrational modes appear. The three normal modes essentially correspond to 
relative center of mass motion with the lowest energy, isovector and isoscalar motion 
with out-of-phase and in-phase oscillations of the two subsystems as the second and 
third mode, respectively. Only the highest-lying isoscalar vibration is present in one- 
component systems, where the excitation energy of -\/5 times the energy quantum of 
the harmonic trap almost precisely is found in our computations. 

New stability conditions are established. The relative center of mass degree of 
freedom provides both more stable structures and reduced stability compared to mean- 
field calculations. Even when the independent subsystems are marginally stable, a 
comparably strong repulsion between the subsystems destabilizes the total system. 

In conclusion, new features arise for two-component systems, i.e. different stable 
structures, excitation modes, and stability conditions. 

Appendix A. The interaction Vab 

We want to compute the expectation value of the interaction Vab in equation (0) over 
all angles, i.e. all coordinates except pA, Pb and r. We assume that the constant angular 
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wave function $o is normalized in the corresponding space, i.e. J (ir2|$oP = 1) where 
dfl = dflrdfl]\f^^i J dQjyg^i with the notation from equation (fTT|) precisely defined in 
[20] • The boson symmetry implies that this expectation value is 

^ Na N 

{%\Vab\%) = 9AB dQj2 Yl S'^'\r,-fj)\%\' 

i=l j=NA+l 

= gABNANB\M^ j dn6^'\fr,^ - rWJ , (A.l) 

where the position vectors can be expressed by center of mass coordinates and Jacobi 
vectors ff, i.e. 



= Ra + \l — (A.2) 



Nb - 1 ^ , , , 

= Rb + ^ VNb-1 ■ (A-3) 

With tat^ — r^B = ^ri — r, we need the 5-function in spherical coordinates, i.e. 

- f) = ^(!Vzl) y ?^±lp^(cos^), (A.4) 

rr/T ^ 4:71 

where P£(cos 6) are Legendre polynomials, 6 is the angle between and f, and = |fVj|. 
We then get 



{%\Vab\%) = QabNaNbI^I^ J dn^A-idnNs- 



S{rn - r) 



(A.5) 



where can be rewritten, by using rj^A-i = Pna-i^^^(^Na-i VNb-i 
PNb-i sinaNg-i, as 

r = —- p^sm aiv^_i + — sm onb-i 

i\a i^b 



- 2W d paPb smaTv^.i smaAr^„i cosfc^Ai? (A. 6) 

where 6ab is the angle between fj^j^_i and tinb-i- 

By changing integration variable from 6ab to a; = — r we arrive at equation (|26|) 
where the reduced interaction is defined by 

I mANBipAPB)'^' / 2r(^)r(M|^) xi/2 
r(iv^-i)(iVs-i)Ur(^)r(3^)/ ^ 

in terms of the basic integral 
I{pA,PB,r) = J daNA-i^^^o^NA-icos^^'^''^ ona-i 
X y (iaAr^_i sinajv^.i cos^^^~'''aiVg_i 

xe(rj+) -r)e(r -rj")) (A.8) 



Two- component boson systems with hyperspherical coordinates 
where © is the Heaviside function and 



-pAsinaNA-i ± 



1 



Pb sm aNg-i] 
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(A.9) 



Na ^'A . y 

fti |pAsinaAr^_i ± Pb sin aNB-i\ , 

where we from now on shall use the last approximation valid for ^ 1 and Nb ^ 1. 

For further derivation we treat separately various cases of different relative sizes 
Pa, Pb and r. The resulting closed analytic expression are one- dimensional integrals of 
simple functions of the angle a, i.e. 



f±{PA,PB,r, a) 



sm a 



cos' 



.3Na-7 



a[l 



±PA sina)^ ^(3jv^_6y2 



3Nb-6 

The integration limits are defined by the functions 
9±{pA: Pb: r) = arcsin[|r ± Pb\I Pa] ■ 

We distinguish between a number of cases. 

1. [2pB < Pa and < r < ps] or [pB < Pa < '^Pb and Q <r < pA- Pb\- 

r9+ rg- 
IipA,PB,r) = / /- - / /+ ■ 
Jo Jo 

2. 2pB < Pa and Pb < r < pa - Pb- 

r9+ 

I(pA,PB,r) = / /_ . 

3. 2pB < Pa and Pa- Pb <r < pa + Pb- 

i-k/2 



I(pA,PB,r) = 



4. 2pB < Pa and Pa + Pb < r: 

I(pA,PB,r) = . 

5. Pb < Pa< '^Pb and Pa - Pb < r < pb- 

I{pA,PB,r) = /- - / /+ ■ 

JO ^0 

6. Pb < Pa< '^Pb and Pb <r < Pa + Pb- 

I{pA,pB,r) = /_ . 

The cases not covered are obtained by interchanging pA and pb- 
The limit of r = is obtained by expansion. We first assume pa > Pb 



r 

dl 

dr 



dr 

2pA 



+ 



OH 



r=0 



r=0 



PB Jo 



6 dr^ 

arcsin(pB/pA) 



r=0 



da sin^ a 



X cos' 



Pb 



(A.IO) 
(A.ll) 

(A.12) 
(A.13) 

(A.14) 
(A.15) 
(A.16) 

(A.17) 



(A.18) 



(A.19) 
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;>arcsin{pfl/pA) 

sin^ a 

pR Jo 



r=0 



Pb Jo 

X cos^^^"^ a[l - 4 sin' a]^^^^^-^")/' (A.20) 
Pb 

3 /•arcsm{pB/pA) 

+ 2(3A^B-8)(3A^ij- 10)^ / c/asin^a 



^2 



where 



r=0 



X cos=^^^-^a[l - ^sin2a](3^^-i2)/2^ 21) 



0. The expressions for > p^ are found by interchanging. 



Appendix B. Stability of p-modes 

When we assume ttia = uib = m, Na = Nb, Pa = Pb and r = the eigenvalues for the 
matrix in equation (|42j) are simply ± The lowest eigenvalue is negative 

when 

(7 p^ = ^ ^ + ^ (aA - ^ < B.l 

4 mp^ TT^/"^ mpj^ 4 

for the solution p^ = p^min obtained from equation (jSTj). This is never possible when 
(^A > O'Ab/'^ and the then the system is always stable. When a a < clab we continue by 
subtracting the function /(pAmin) = in equation (jSZ|) from the inequality in equation 
(jB.ljl . This immediately leads to 

(3iV2)V2 5 • 3V2iv3/2 _ 

^Y7^ < PAmin < 227rl/2 = ' ^^'"^^ 

where we added the inequality for the lowest possible value of pAmin- Then g{pAb) < 
is the condition for instability, which by simple insertion results in equation PHj) . Using 
instead the lower limit in equation ()B.2|) we arrive at condition which is the stability 
condition for a one-component system. 



References 

[1] Pethick C J and Smith H 2001 Bose-Einstein Condensation in Dilute Gases (Cambridge: 
Cambridge University Press) 

[2] Pitaevskii L and Stringari S 2003 Bose-Einstein Condensation (Oxford: Clarendon Press) 

[3] Myatt C 3 et al 1997 Phys. Rev. Lett. 78 586 

[4] Hall T)S et al 1998 Phys. Rev. Lett. 81 1539 

[5] Maddaloni P et al 2000 Phys. Rev. Lett. 85 2413 

[6] Modugno G et al 2002 Phys. Rev. Lett. 89 190404 

[7] Truscott AG et al 2001 Science 291 2570 

[8] Schreck F et al 2001 Phys. Rev. Lett. 87 080403 

[9] Roati G et al 2002 Phys. Rev. Lett. 89 150403 

[10] Modugno G et al 2002 Science 297 2240 

[11] Hadzibabic Z et al 2002 Phys. Rev. Lett. 88 160401 

[12] Chui S T and Ao P 1999 Phys. Rev. A 59 1473 



Two- component boson systems with hyperspherical coordinates 



Esry B D and Greene C H 1999 Phys. Rev. A 59 1457 
Gordon D and Savage C M 1998 Phys. Rev. A 58 1440 
Pu H and Bigelow N P 1998 Phys. Rev. Lett. 80 1134 
Ohberg P 1999 Phys. Rev. A 59 634 

Svidzinsky A A and Chui S T 2003 Phys. Rev. A 68 013612 
Busch Th et al 1997 Phys. Rev. A 56 2978 
Bohn et al 1998 Phys. Rev. A 58 584 
S0rensen O et al 2002 Phys. Rev. A 66 032507 
S0rensen O et al 2002 Phys. Rev. Lett. 89 173002 

Sogo T et al 2004 submitted for publication [Preprint cond-mat / 040227l| | 
S0rensen O et al 2004 J. Phys. B 37 93 
S0rensen O et al 2003 Phys. Rev. A 68 063618 



